function knearestneighbour_mean_interpolator, gridme, gridsize, k, stages=stages
	bounds = findbounds(gridme)
	return_me = dblarr(gridsize, gridsize)
	coefficients = rebin(dindgen(gridsize) / gridsize, gridsize, gridsize)
	cellx = bounds[0] + (bounds[1] - bounds[0]) * coefficients
	celly = bounds[2] + (bounds[3] - bounds[2]) * transpose(coefficients)
	cellcentres = transpose([[cellx[*]], [celly[*]]])
	cellsz = size(cellcentres, /dimensions)
	;for big data sets must do in several stages
	stage = (keyword_set(stages) ? stages : 1)
	icnt = 0l
	for i=0, stage - 1 do begin
		print, "Begining stage ", i
		kindices = knearestneighbour(cellcentres[*,i*(cellsz[1]/stage):i*(cellsz[1]/stage) + (cellsz[1]/stage) - 1], gridme[0:1,*], k)
		for j=0, n_elements(kindices[0,*]) - 1 do begin
				p =  array_indices(return_me, icnt++)
				return_me[p[0],p[1]] = mean(gridme[2,kindices[*,j]])
		endfor
		;for j=0, cellsz[1]/stage - 1 do $
		;	return_me[(cellcentres[0,i:i+(cellsz[1]/stage)-1])[j], (cellcentres[1,i:i+(cellsz[1]/stage)-1])[j]] = mean(gridme[2,kindices[*,j]])
	endfor
	return, return_me
end

function knearestneighbour_mean_interpolator_old_AND_COOL, gridme, gridsize, k, stages=stages
	bounds = findbounds(gridme)
	return_me = dblarr(gridsize, gridsize)
	coefficients = rebin(dindgen(gridsize) / gridsize, gridsize, gridsize)
	cellx = bounds[0] + (bounds[1] - bounds[0]) * coefficients
	celly = bounds[2] + (bounds[3] - bounds[2]) * transpose(coefficients)
	cellcentres = transpose([[cellx[*]], [celly[*]]])
	cellsz = size(cellcentres, /dimensions)
	;for big data sets must do in several stages
	stage = (keyword_set(stages) ? stages : 1)
	for i=0, stage - 1 do begin
		kindices = knearestneighbour(cellcentres[*,i:i+(cellsz[1]/stage)-1], gridme[0:1,*], k)
		for j=0, gridsize - 1 do $
			for l=0, gridsize - 1 do $
				return_me[l,j] = mean(gridme[2,kindices[*,l*j]])
		;for j=0, cellsz[1]/stage - 1 do $
		;	return_me[(cellcentres[0,i:i+(cellsz[1]/stage)-1])[j], (cellcentres[1,i:i+(cellsz[1]/stage)-1])[j]] = mean(gridme[2,kindices[*,j]])
	endfor
	return, return_me
end